Q 

AD-A218  144 

ICASE 

TRIANGLE  BASED  TVD  SCHEMES  FOR  HYPERBOLIC 
CONSERVATION  LAWS 


flflC  FILE  COZY 

NASA  Contractor  Report  181984 
ICASE  Report  No.  90-10 


Louis  J.  Durlofsky 
Stanley  Osher 
Bjorn  Engquist 


Contract  No.  NAS  1-18605 
January  1990 


Institute  for  Computer  Applications  in  Science  and  Engineering 
NASA  Langley  Research  Center 
Hampton,  Virginia  23665-5225 

Operated  by  the  Universities  Space  Research  Association 


NASA 

National  Aeronautics  and 
Space  Administration 

Langley  Research  Center 

Hampton,  Virginia  23665-5225 


[ 


DanuBtmoN  statement  a 

Approved  for  public  rslsosei 
_ Distribution  Uaiiautsd 


90  02  21  006 


TRIANGLE  BASED  TVD  SCHEMES  FOR  HYPERBOLIC 

CONSERVATION  LAWS 


Louis  J.  Durlofsky 

Chevron  Oil  Field  Research  Company 
P.O.  Box  446 
La  Habra,  CA  90633-0446 

Stanley  Osher1’2  and  Bjorn  Engquist2 
Department  of  Mathematics 
UCLA 

Los  Angeles,  CA  90024 

ABSTRACT 


A  triangle  based  TVD  (total  variation  diminishing)  scheme  for  the  numerical  approxima¬ 
tion  of  hyperbolic  conservation  laws  in  two  s--ace  dimensions  is  constructed.  The  novelty  &f 
the  scheme  lies  in  the  nature  of  the  preprocessing  of  the  cell  averaged  data,  which  is  accom¬ 
plished  via  a  nearest  neighbor  linear  interpolation  followed  by  a  slope  limiting  procedure. 
Two  such  limiting  procedures  are  suggested.  The  resulting  method  is  considerably  more 
simple  than  other  triangle  based  non-oscillatory  approximations  which,  like  this  scheme,  ap¬ 
proximate  the  flux  im  to  second  order  accuracy.  Numerical  results  for  linear  advection  and 
Burgers’  equation  am  presented. 
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1.  Introduction 


In  the  last  ten  years  there  has  been  considerable  effort  aimed  at  constructing 
and  analyzing  high  order  accurate,  non-oscillatory  approximations  to  hyperbolic 
conservation  laws  (see  e.g.,  [2],  [8]).  It  is  by  now  well  established  that  the  spon¬ 
taneous  development  of  shock  waves  and  the  appearance  of  steep  gradients  in  the 
solution  require  higher  order  schemes  to  have  an  adaptive  stencil  (by  adaptive  stencil 
we  mean  an  adaptive  flux  approximation,  not  an  adaptive  grid)  in  order  to  suppress 
the  spurious  oscillations  that  plague  conventional  finite  difference  methods.  Total 
variation  diminishing  (TVD)  schemes,  one  such  class  of  second  order  accurate  meth¬ 
ods  that  eliminate  unphysical  oscillations,  have  been  used  successfully  in  a  variety 
of  applications.  Recently,  a  new  class  of  methods,  essentially  non-oscillatory  (ENO) 
schemes  ([3],  [4]),  which  surpass  the  second  order  accurate  barrier  associated  with 
TVD  schemes,  has  been  developed.  An  alternative  approach  for  third  order  schemes 
was  developed  in  [10]. 

Extensions  of  TVD  and  ENO  schemes  to  two  and  three  dimensions  are  typically 
accomplished  in  a  dimension  by  dimension  fashion,  via  space-operator  splitting. 
Therefore,  the  extension  of  these  higher  order  schemes  to  the  solution  of  hyperbolic 
conservation  laws  on  unstructured  grids,  such  as  a  triangular  mesh,  is  not  imme¬ 
diate.  It  is  our  intent  in  this  paper  to  devise  a  second  order  accurate  scheme  of 
TVD  type  which  is  applicable  to  an  unstructured  triangular  grid.  Our  scheme  is 
based  on  a  finite  volume  type  discretization  and  is  particularly  straightforward  to 
implement.  The  scheme  relies  on  a  very  local  adaptive  interpolation  idea,  which 
results  in  computational  efficiency.  In  the  future,  we  expect  to  extend  the  adaptive 
two  dimensional  interpolation  ideas  presented  here  to  develop  triangle  based,  higher 
order  ENO  schemes. 

Several  approaches  for  the  solution  of  hyperbolic  conservation  laws  on  trian¬ 
gular  grids  already  exist.  These  techniques  are,  however,  in  the  context  of  finite 
element  methods  and  have  utilized  flux  corrected  transport  (FCT)  ideas  [6]  or  have 
required  the  generation  of  a  complex  auxiliary  grid  [9]  or  are  truly  finite  element 
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methods  in  space  and  time  and  thus  are  more  costly  computationally  [5].  The 
methodology  presented  here  is,  in  our  opinion,  simpler  and  more  efficient,  primarily 
because  a  finite  volume  rather  than  a  finite  element  approach  is  used,  thus  avoiding 
the  overhead  associated  with  finite  element  schemes. 

The  TVD,  second  order  accurate  methods  we  shall  develop  in  §2  are  technically 
neither  total  variation  diminishing  nor  strictly  second  order  accurate.  We  follow  the 
convention  of  calling  two  dimensional  schemes  TVD  if  they  are  formal  extensions 
of  one  dimensional  TVD  schemes,  as  our  scheme  is.  In  general,  however,  the  total 
variation  may  increase  [1],  though  a  maximum  principle  is  satisfied.  Also,  although 
the  fluxes  are  approximated  up  to  second  order,  the  truncation  error  is  technically 
lower  because  of  the  adaptive  stencil  and  the  variable  size  of  the  triangles.  Numerical 
experiments  presented  in  §3  indicate  orders  of  accuracy  between  1.6  and  1.9  in  the 
L\  norm.  In  §4  we  suggest  further  extensions  of  the  method  within  the  TVD  context 
and  indicate  partial  extensions  to  include  diffusive  terms. 

2.  Construction  of  the  Numerical  Schemes 

Our  intent  in  this  section  is  to  develop  a  scheme  to  solve  hyperbolic  conserva¬ 
tion  laws  on  triangular  grids  in  two  space  dimensions.  The  method  presented  is  for 
single  hyperbolic  conservation  laws,  though  hyperbolic  systems  can  be  treated  anal¬ 
ogously  in  a  field  by  field  manner.  Our  method  is  finite  volume  based  and  achieves 
greater  than  first  order  accuracy  through  use  of  a  novel  adaptive  flux  interpolation 
procedure.  We  first  present  the  general  finite  volume  approach,  then  introduce  our 
general  limiting  procedure,  and  then  discuss  various  specific  limiters. 

2.1  Finite  Volume  Discretization 

Consider  the  hyperbolic  conservation  law, 

ut  +  V  •  F(u)  =  g(x,t), 
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u(x,0)  =  u0(x), 


(2.1) 


subject  to  boundary  conditions.  We  wish  to  solve  (2.1)  on  a  triangular  grid,  a 
portion  of  which  is  shown  schematically  in  Fig.  1.  Integrating  (2.1)  over  a  triangle 
(Aabc  to  he  specific)  gives, 


-  / 

dt  L.„c 


udA 


-L 


(V  •  F)dA, 


A^flc 


(2.2) 


where  A abc  represents  both  the  region  ABC  and  its  area,  and  <z(x,t)  has  been 
taken  to  be  zero  for  simplicity  of  exposition  only.  Applying  the  divergence  theorem 
to  the  right  hand  side  of  (2.2)  and  defining, 


u  =  {  udA)(&ABC )  1 , 

J  AytflC 

.e.,  u  is  the  average  of  u  over  A  abc,  gives, 

— [/  F  n  abM+  [  FtiAcdC, 
&ABC  [J(AB  JlAC 

+  /  F  •  rifled^  • 

JeBc  J 


d_ 

dt 


(2.3) 


Note  that  u  is  equal  to  the  value  of  u  evaluated  at  the  triangle  centroid  ( Xabc ) 
to  within  O(Aabc),  or,  analogously,  to  within  0((2),  where  £  is  the  characteristic 
length  of  a  side  of  A  abc-  Here  n  is  the  unit  outward  normal. 


We  approximate  (2.3)  by  first  using  a  semi-discrete  approach  where  the  ap 
proximation  is 

VABc(t)  UABc(t)\ 

the  same  is  true  for  all  triangles.  First  order  accurate  monotone  schemes  can  eas¬ 
ily  be  constructed  -  see  e.g.,  [7],  [14].  Let  hsc(w\  ,w2)  he  a  twr  point  Lipschitz 
continuous  monotone  flux,  approximating  F  •  n bc,  he., 

(2.4a)  hgc(iv,w)  =  F  •  n bc, 

(2.4b)  hgc(wi ,  W2)  is  a  nondecreasing  function  of  iv\  and  a  nonincreasing  func¬ 
tion  of  U>2- 
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Then  our  semidiscrete  monotone  approximation  is, 

h.Bc(vABC,VBCD )  '  ?BC 
+  hAB(vABC,VABF)  •  ?AB  (2-5) 

+  hAc(vABC,  VACE)  '  ?AC  , 

where  £bc  is  the  length  of  the  side  BC ,  etc.  “E”  schemes  may  also  be  used  -  see 
[7]  for  the  definition  and  for  examples. 

To  obtain  higher  order  accuracy  we  preprocess  our  initial  data  so  that  in  each 
triangle,  in  particular  Aabc ,  a  linear  function  is  obtained  whose  cell  average  equals 
vABCi  but  which  is  within  0(A)  of  uabc  in  regions  of  smoothness.  Here  A  is  the 
maximum  area  of  the  four  triangles  seen  in  Fig.  1.  Moreover,  this  linear  function 
will  not  introduce  new  oscillations  in  our  approximation.  This  (simple)  construction 
is  the  key  part  of  this  paper;  it  will  be  described  at  the  end  of  this  section.  We  call 
this  linear  approximation  L&(x).  It  is  generally  discontinuous  across  the  boundary 
of  each  triangle. 

Let  xbc  be  the  midpoint  of  side  BC ,  etc.  Let  La(  XBc)  denote  the  limit  of 
La(x)  as  x  — »  Xbc  from  inside  triangle  ABC  and  L&(x0BC)  denote  the  limit  as 
x  — »  xbc  from  outside  triangle  ABC.  Generally, 

I  La(x1Dc)  —  £a(xbc)I  =  0(A). 

Our  TVD,  second  order  accurate,  semi-discrete  approximation  to  (2.3)  is 

a  ,  x  i 
aiVABc(t)  =  “  w 

+  hAB{L&(xlAB),  L<\(x°ab))  •  £ab  (2-6) 

+  hAc{L^{x\c),  La.(x°ac))  •  l AC  • 

By  the  midpoint  formula  for  integrals,  this  approximation  is  weakly  second  order 
accurate,  in  the  sense  that  each  of  the  three  flux  terms  above  is  within  0(A)  of 


h-Bc{L^{xlBC),  La(x°bc))  ■  £bc 


^VABc(t)  =  - 
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the  line  integrals,  f  F  ■  ruff,  along  the  corresponding  interfaces.  However  due  to 
the  shifting  stencil  and  varying  size  and  relation  of  the  triangles,  the  pointwise 
truncation  error  is  generally  only  0(A  2 ),  i.e.,  first  order.  The  performance  appears 
to  be  around  1.6-1. 9  order  in  L\  for  smooth  flow  (see  §3). 

2.2  Construction  of  Linear  Function  La 

We  now  describe  the  construction  of  La-  In  each  interior  triangle,  three  can¬ 
didates  for  La ,  designated  L\ ,  are  generated.  The  first  such  candidate  is  the 
linear  interpolate  of  the  three  values 

(xabc,vabc),  ( xbcd,vbcd ),  ( *-ace,vace ), 

Ll  is  the  interpolation  of 

( xabc,vABc ),  ( *bcd,vbcd ),  ( *abf,vabf ), 

and  L\  the  interpolation  of 

(x.ABC,VABc),  (*-ace,vace),  (x-abf,vabf)- 

These  three  linear  interpolants  are  sketched  in  Fig.  2.  Here  and  below  we  assume 
that  the  three  triangle  centroids,  xabc,*bcd  and  X-abf  are  not  colinear.  At  this 
point,  three  possible  L ^  exist,  and  a  limited  version  of  La  must  be  selected  from 
these.  To  accomplish  this,  we  first  compute  the  magnitude  of  the  gradient  of  each 
L'a i  i.e., 

A  A 

[(— L\)2 +  {—L\f}'  =|VZ4|,  for  i  =  1,2,3.  (2.7) 

By  analogy  with  limiting  procedures  in  one  space  dimension  ([13]),  a  valid,  though 
very  non-compressive  limiter,  corresponds  to  the  selection  of  the  L\  for  which  |VL^| 
is  the  minimum.  This  choice  is  analogous  to  the  min  limiter  in  second  order  ENO 
methods  ([3]);  no  special  precautions  need  be  taken  at  extrema. 
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It  is  desirable  to  construct  a  more  compressive  limiter  than  that  described 
above,  particularly  for  problems  involving  linear  or  contact  discontinuities.  To  ac¬ 
complish  this,  we  first  consider  the  more  compressive  slope  limiters  in  one  dimension, 
the  $  type  limiters  described  by  Sweby  [13]  (his  equation  3.17)  of  which  superbee 
is  the  most  compressive,  corresponding  to  $  =  2.  These  limiters  allow  the  use  of 
piecewise  linear  approximations  to  the  solution  for  which  the  slope  is  not  the  min¬ 
imum,  subject  to  the  restriction  that  no  overshoot  (or  undershoot)  occurs  at  the 
cell  boundaries. 

The  next  limiter  we  describe  is  a  multidimension  analog  of  the  one  dimensional 
$  limiters.  The  approach  here  is  to  select  the  L'A  for  which  |V£^|  is  maximized, 
subject  to  the  restriction  that  no  overshoot  or  undershoot  occurs  at  any  of  the  three 
triangle  boundaries.  The  procedure  is  as  follows: 

(i)  Select  the  L'A  for  which  |V£}J  is  the  maximum. 

(ii)  Check  for  overshoot  or  undershoot  at  Xab,  *AC  and  Xbc •  For  L'A  to  represent 
a  valid  £a,  it  suffices  to  verify  that,  for  A abc, 

L&{xac)  is  between  vA bc  and  vace, 

L&(xab)  is  between  vabc  and  vabf  and 

Ta(xbc)  is  between  vAbc  and  vbcd- 

If  these  three  requirements  are  satisfied,  L'A  is  the  appropriate  L&. 

(iii)  If  the  L\  above  results  in  overshoot  or  undershoot  at  any  one  of  the  three 
midpoints,  select  the  LlA  for  which  ]V£^j  is  the  second  largest  and  repeat  the 
test  in  (ii).  If  this  L\  does  not  satisfy  the  test  in  (ii),  select  the  L\  for  which 
|V£]J  is  the  minimum  and  again  proceed  through  the  test  in  (ii). 

(iv)  If  all  L\  fail  (ii),  revert  to  a  piecewise  constant  approximation  for  A  abc',  he., 
L  a  =  vabc- 

Given  £a,  the  right  hand  side  of  (2.6)  can  be  evaluated  and  vABc(i )  integrated 
in  time.  This  time  integration  is  accomplished  via  a  second  order  TVD  Runge-Kutta 
procedure  [11]. 


6 


3.  Numerical  Verification  of  Higher  Order  Scheme 


In  this  section  we  present  results  for  the  convergence  of  the  general  method 
described  in  §2,  as  well  as  solution  contours  and  profiles  demonstrating  the  accuracy 
of  the  method.  In  all  cases,  the  solution  region  is  a  square  domain  discretized  via 
right  triangular  ‘volumes’  (referred  to  as  elements),  as  shown  in  Fig.  3.  Periodic 
boundary  conditions  are  imposed  in  both  the  x—  and  the  y—  directions;  the  initial 
condition  is  similarly  x—  and  y  —  periodic.  In  all  cases  the  more  compressive  limiter 
described  in  §2.2  is  used  for  the  TVD  scheme. 

3.1  Rate  of  Convergence 

To  assess  rate  of  convergence,  the  scheme  is  applied  to 
conservation  law 

ut  +  V  •  (au)  =  0, 

subject  to  the  initial  condition 

u0(x,y)  =  sin(27rx)sin(27ry).  (3-2) 

Our  base  monotone  scheme  uses  the  EO  flux  [7]: 


h(wuw2)  =  /+(u>i)  +  f-(w2).  (3.3) 

For  linear  equations  with  constant  a  =  (ax,ay), 

f+{u)  =  [max((a  •  n),0)]u,  (3.4a) 

f-(u)  =  [min((a  ■  n),0)]u.  (3.46) 

For  Burgers’  equation  (considered  below),  where  f\  =  f2  =  (1/2 )ti2, 

f+(u)  =  max((n*  +  Uy  )u2, 0),  (3.5a) 

/_(u)  -  min((rc*  y  Hy  )tt2, 0),  (3.56) 
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the  solution  of  the  linear 
(3.1) 


where  nx  and  ny  represent  the  components  of  n. 


A  contour  plot  of  the  initial  condition  (3.2)  is  shown  in  Fig.  4.  Four  extrema 
are  evident.  The  rate  of  convergence  of  the  method  was  determined  for  both  the 
case  ax  =  ay  =  1  and  ax  =  1,  ay  =  0.  Further,  convergence  was  assessed  both  on 
an  element  by  element  basis  and  after  applying  a  local  averaging  procedure.  It  is 
expected  that  local  averaging  procedures  would  enhance  the  rate  of  convergence,  as 
the  scheme  is  expected  to  be  second  order  in  only  the  weak  sense;  i.e.,  after  integrat¬ 
ing  locally  in  space  and  time  (a  type  of  local  averaging).  The  averaging  performed 
in  this  study  is,  however,  only  spatial;  no  temporal  averaging  is  attempted.  This 
is  because  spatial-temporal  averages  are  rather  cumbersome  to  perform  in  practice, 
and  the  spatial  averaging  alone  reveals  the  expected  trend.  Computations  were 
performed  for  grids  ranging  in  discretization  from  200  elements  (l  —  0.1,  where  £ 
is  the  spacing  between  adjacent  nodes  or,  analogously,  i  —  (2A)1/2,  with  A  the 
area  of  any  element)  to  12800  elements  (£  =  0.0125).  In  all  cases  the  CFL  number, 
A(=  A t/£),  was  set  to  0.1. 

Displayed  in  Fig.  5  is  a  log-log  plot  of  L\  error  versus  i.  In  this  case,  az  =  ay  = 
1.  Results  are  shown  for  both  a  first  order  scheme  and  the  higher  order  scheme, 
with  error  computed  on  an  element  by  element  basis.  Least  squares  linear  fits  give 
the  order  of  convergence  for  the  two  methods;  for  the  first  order  method  we  obtain 
0.93  and  for  the  higher  order  method  1.77.  Figure  6  displays  an  analogous  plot 
after  applying  a  local  averaging  procedure.  Specifically,  this  averaging  procedure 
entails  averaging  the  computed  value  of  u  over  square  regions  comprised  of  two 
adjacent  elements  and  computing  the  error  in  terms  of  the  difference  between  this 
average  and  the  exact  solution  of  Eq.  (3.1)  evaluated  at  the  square  midpoint.  For 
the  grid  displayed  in  Fig.  3,  100  such  square  regions  exist.  Again,  averaging  is 
only  applied  spatially;  no  temporal  averaging  is  performed.  Assessing  error  in  this 
manner  results  in  least  squares  linear  fits  of  slope  0.94  for  the  first  order  method  and 
1.81  for  the  higher  order  method.  As  expected,  local  averaging  enhances  the  rate 
of  convergence  though,  in  this  case,  the  improvement  is  minimal.  In  other  cases, 
however,  the  improvement  is  more  substantial.  For  example,  using  ar  =  1,  ay  =  0 
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in  Eq.  (3.1)  yields  the  following  convergence  results.  For  the  first  order  method, 
convergence  is  0(0.99)  with  no  local  averaging  and  0(1.06)  with  averaging.  For  the 
higher  order  scheme,  the  convergence  rates  are  0(1.60)  and  0(1.73),  respectively. 

Results  for  L2  and  Loo  error  display  slower  rates  of  convergence.  For  the  case 
ax  =  ay  —  1,  L2  error  is  ~  0(f164)  with  local  averaging  and  0(£162)  on  an  element 
by  element  basis  with  Loo  error  ~  0(^°-94)  with  local  averaging  and  0(£°'91)  with  no 
averaging.  The  expected  result  for  L2  error  is  ~  0(E15)  and  for  L0 0  error  ~  0(E),  as 
in  one  dimensional  TVD  methods.  Although  the  discrepancies  between  the  expected 
and  numerical  results  are  relatively  slight,  it  is  not  clear  why  the  L2  error  converges 
faster  than  expected  while  the  L0 0  error  converges  slower  than  expected. 

Shown  in  Table  1  is  a  compilation  of  the  rates  of  convergence  of  L\ ,  L2  and  Loo 
error.  Results  for  both  a  first  order  scheme  and  our  more  compressive  TVD  scheme 
are  displayed.  In  all  cases  the  initial  condition  is  as  in  (3.2).  Error  is  computed 
over  the  entire  domain  in  two  ways:  (1)  element  by  element  and  (2)  by  combining 
two  adjacent  elements  into  squares.  In  all  cases  E  ranges  from  0.0125  to  0.1. 

Slightly  improved  rates  of  convergence  in  L\  are  obtained  when  the  initial 
condition  contains  no  extrema.  This  is  demonstrated  in  Table  2,  where  results 
for  L\  error  for  the  initial  condition  uo(x,y )  =  sin(7rx/2)  sin(7rj//2)  are  displayed. 
Here,  to  eliminate  the  effects  of  the  discontinuity  in  u  at  the  boundary  (recall  that 
periodic  boundary  conditions  are  imposed),  error  is  computed  only  over  the  region 
0-6  <  x,y  <  0.8  at  an  early  time,  t  =  0.05.  In  one  case,  local  averaging  has  a  more 
dramatic  effect,  improving  the  L\  accuracy  of  the  TVD  scheme  from  1.22  to  1.80. 

Based  on  the  numerical  results  presented  above  and  the  analysis  presented  in 
§2,  we  feel  that  the  method  can  be  considered  to  be  second  order  accurate  in  L\  in 
the  weak  sense.  Though  our  convergence  results  always  indicate  convergence  slower 
than  quadratic,  this  is,  in  our  opinion,  due  to  the  fact  that  these  results  are  not 
strictly  measuring  weak  convergence.  If  such  convergence  could  be  unambiguously 
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measured,  it  is  our  contention  that  the  method  would  indeed  display  second  order 
convergence. 

3.2  Examples  of  Numerical  Accuracy 

We  now  present  some  detailed  numerical  results  for  our  second  order  scheme 
and  compare  these  with  the  results  of  a  first  order  method.  The  first  results  are  for 
the  solution  of  Eqs.  (3.i)  and  (3.2)  with  az  =  ay  =  1.  Figure  7  displays  the  solution 
contour  results  for  the  first  order  scheme  with  800  elements  (£  =  0.05)  and  A  =  0.1 
(the  same  CFL  number  is  used  in  all  computations)  at  t  =  1.  The  exact  solution 
is  a  reproduction  of  the  initial  condition,  shown  in  Fig.  4.  The  first  order  method 
is  clearly  very  diffusive;  the  maximum  value  of  u  is  here  only  0.25,  in  contrast  to 
the  maximum  in  the  initial  condition  of  1.  Results  for  the  second  order  scheme 
at  t  =  1  are  shown  i  l  Fig.  8.  Though  some  distortion  of  the  initial  condition  is 
apparent,  the  solution  is  considerably  improved  over  the  first  order  solution;  the 
maximum  value  of  u  is  now  0.76.  Shown  in  Fig.  9  are  the  t  —  1  results  for  the 
first  order  scheme  using  3200  elements  (£  =  0.025).  Substantial  numerical  diffusion 
is  still  evident;  the  maximum  value  of  u  is  only  0.49.  The  solution  contour  using 
the  second  order  method  is  displayed  in  Fig.  10.  The  t  —  1  solution  in  this  case 
closely  resembles  the  initial  condition,  with  a  maximum  value  of  u  of  0.88.  Figures 
11  and  12  show  solution  profiles  taken  along  the  line  y  =  x  (the  velocity  direction) 
at  t  =  0,  0.25,  0.5,  0.75  and  1  for  both  the  first  and  second  order  methods.  In 

both  cases,  3200  elements  were  used.  The  second  order  results  are  quite  sharp  at 
all  times,  while  the  first  order  results  show  a  continual  degradation  with  increasing 
time. 

Solution  profiles  for  computations  using  the  second  order  method  with  a  non¬ 
linear  flux  function,  /  =  (l/2)u2  in  Eq.  (2.1)  (he.,  the  inviscid  Burgers’  equation), 
with  the  initial  condition  u0(x ,y)  =  sin(27rx),  are  shown  in  Fig.  13.  Though  this 
is  an  essentially  one  dimensional  problem,  no  overshoot  or  unphysical  oscillations 
appear  in  the  solution.  The  solution  is  sharper  with  the  second  order  method  than 
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with  a  first  order  scheme  (first  order  results  are  not  shown),  though  the  improvement 
is  of  course  less  dramatic  than  in  the  linear  examples  presented  above. 

4.  Possible  Extensions 

Other  limiting  procedures  are  quite  feasible  and  should  be  tested.  Our  com¬ 
pressive  limiter  is  not  a  direct  analogue  of  superbee,  since  superbee  (and  many  other 
limiters  [13])  occasionally  allows  values  other  than  zero  or  any  of  the  slopes  being 
compared  to  be  the  final  choice  of  slope  (or  gradient  in  our  two  dimensional  case). 

A  more  significant  issue  is  the  treatment  of  diffusive  terms.  In  this  case,  the 
governing  equation  is  of  the  form 

ut  +  V  ■  F(tt)  =  e(uxx  +  uyy),  e  >  0.  (4.1) 

The  discrete  analogue  of  (2.3)  now  involves  the  additional  term, 


on  the  right  side  of  (2.3).  Up  to  first  order  accuracy,  we  compute  each  of  the  three 
terms  in  (4.2)  as  follows.  The  limiting  procedure  has  already  given  us  a  gradient 
within  the  triangle  ABC  as  well  as  for  each  of  the  three  neighbors.  Therefore,  the 
integral  along  side  AB  in  (4.2)  can  be  computed  approximately  as 

l(VLABC  +  VLABF)-n ]^£.  (4.3) 

The  integrals  along  the  other  sides  are  approximated  analogously.  This  is  generally 
a  first  order  accurate  method  (second  order  accuracy  occurs  in  special  cases;  e.g.,  if 
all  the  triangles  are  equilateral).  However,  since  e  is  relatively  small  here  (otherwise 
transport  is  diffusion  dominated  and  the  sophisticated  treatment  of  convection  is 
unnecessary),  we  believe  this  to  be  an  adequate  treatment  of  these  terms. 

Finally,  we  mention  that  '.vork  is  underway  to  approximate  (2.1)  using  a  higher 
order  accurate  ENO  triangle  based  method.  See  [11],  [12]  for  successful  Cartesian 
coordinate  approaches. 
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Table  1.  Computed  accuracy  of  TVD  scheme  for  the  linear  case.  Initial  condition 
u0(x,y )  =  sin(27ri)  sin(27ry).  Error  computed  over  the  entire  domain. 


az  =  a 

*  =  1 

Scheme 

Norm 

#  elements 

n 

2nd  0 

Ii 

1 

1.77 

2nd  0 

Li 

2 

1.81 

2nd  0 

l2 

1 

1.62 

2nd  0 

l2 

2 

1.64 

2nd  0 

Lqo 

1 

0.91 

2nd  0 

Loo 

2 

0.94 

1st  0 

Lx 

1 

0.93 

1st  0 

Lx 

2 

0.94 

1st  0 

l2 

1 

0.94 

1st  0 

l2 

2 

0.94 

1st  0 

Loo 

1 

0.95 

1st  0 

Loo 

2 

0.94 

ft 

H 

II 

o 

II 

a> 

<3 

Scheme 

Norm 

#  elements 

n 

2nd  0 

Lx 

1 

1.60 

2nd  0 

Lx 

2 

1.73 

1st  0 

Lx 

1 

0.99 

1st  0 

Lx 

2 

1.06 
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Table  2.  Computed  accuracy  ( Li )  of  TVD  scheme  for  the  linear  case.  Initial 
condition  u0(x,y)  =  sin(n/2)  sin(7ry/2)  contains  no  extrema.  Error  computed 
over  0.6  <  x,  y  <  0.8  at  t  =  0.05. 


az  =  ay  =  1 

Scheme  #  elements  n 

2nd  0  1  1.85 

2nd  0  2  1.87 

CX  —  0,  dy  —  1 

Scheme  $  elements  n 

2nd  0  1  1.22 

2nd  0  2  1.80 

1st  O  1  0.99 

1st  0  2  1.10 
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Figure  3 

Triangular  grid  used  for  the  numerical  calculations. 
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Figure  4 

Contour  plot  of  the  initial  condition  (3.2).  Contours  correspond 
to  u  =  0,  ±0.15,  ±0.3,  ±0.45,  ±0.6,  ±0.75,  ±0.9. 
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Figure  5 

L  i  error  on  a  per  element  basis  for  the  case  ax  =  a  y  =  1  for 
first  order  (O)  and  second  order  (x)  schemes.  Lines  are 
least  square  fits  with  slopes  as  indicated. 


20 


X 


Figure  6 

Li  error  after  local  averaging  for  the  case  ax  =  ay  =  1  for 
first  order  (O)  and  second  order  (>4  schemes.  Lines  are 
least  square  fits  with  slopes  as  indicated. 
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Figure  1 1 

Solution  profiles  along  the  line  y=x  for  the  first  order  scheme 
(3200  elements). 
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scheme  (3200  elements). 
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Figure  13 

Solution  profiles  for  Burgers'  equation  using  second  order 
scheme  (800  elements). 
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